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We analyze the structure of stochastic dynamics near either a stable or unstable fixed point, 
where force can be approximated by linearization. We find that a cost function that determines a 
Boltzmann-like stationary distribution can always be defined near it. Such a stationary distribution 
does not need to satisfy the usual detailed balance condition, but might have instead a divergence- free 
probability current. In the linear case the force can be split into two parts, one of which gives detailed 
balance with the diffusive motion, while the other induces cyclic motion on surfaces of constant cost 
function. Using the Jordan transformation for the force matrix, we find an explicit construction 
of the cost function. We discuss singularities of the transformation and their consequences for the 
stationary distribution. This Boltzmann-like distribution may be not unique, and nonlinear effects 
and boundary conditions may change the distribution and induce additional currents even in the 
neighborhood of a fixed point. 
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I. INTRODUCTION 



In equilibrium statistical mechanics an important role 
is played by the principle of detailed balance and by the 
related fluctuation-dissipation theorem. Einstein used 
the principle that the excess energy that is put into each 
mode of an equilibrium system in the course of thermal 
fluctuations is also removed from the same mode by dis- 
sipative forces. This is implicit in his work on Brownian 
movement 0, and explicit in later works on the photo- 
electric effect 0, and on the relation between sponta- 
neous and induced emission of electromagnetic radiation 
0. This principle was formulated as the principle of 
detailed balance by Bridgman and used to explain 
Johnson noise in electrical circuits by Nyquist p|. This 
is related to the fact that the same processes that drive 
fluctuations in the neighborhood of a typical equilibrium 
configuration also drive the configuration back towards 
a typical equilibrium or steady state configuration when 
it is displaced from equilibrium by an amount which is 
small, but large compared with the fluctuations in ther- 
mal equilibrium. In this situation the equilibrium distri- 
bution in phase space is just the Boltzmann distribution, 
proportional to e~@ E , where (3 is inversely proportional 
to temperature, and E is the energy of the point in phase 
space. Configurations that differ significantly from those 
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that contribute to the minimum of the free energy are 
driven back to the neighborhood of this minimum by dis- 
sipative effects such as thermal or electrical conduction, 
or viscosity, and the magnitude of these effects is related 
to the equilibrium fluctuations of related variables. 

In many situations there is no thermodynamic equi- 
librium, but external steady and fluctuating forces drive 
the system into a steady or very slowly varying state for 
which the principle of detailed balance does not hold. A 
light bulb powered by an external battery, or a chemi- 
cal reaction in which the reactants are introduced at a 
steady rate and the products of the reaction removed at 
a steady rate, would both be examples of such a situa- 
tion. Even in a situation which is almost in equilibrium, 
such as a system which is started in equilibrium at a local 
minimum of the free energy, but which can can go over 
a saddle point to a deeper minimum, the behavior near 
the saddle point does not satisfy the principle of detailed 
balance, since there is a current over the saddle point. 

For such systems without detailed balance there is no 
general method of obtaining the equilibrium distribution 
from a knowledge of the steady and stochastic forces, 
such as the Boltzmann distribution provides for a system 
with detailed balance. In recent work one of us [|J has 
developed a method valid near a stable fixed point which, 
even when detailed balance does not hold, obtains a cost 
function analogous to the energy for the Boltzmann dis- 
tribution. If this method can be extended away from the 
linear region in the neighborhood of a fixed point it may 
provide a new method for dealing with problems of this 
sort 0. 

Great efforts have been spent on finding such a cost 



2 



function ever since the work of Onsager ||- Results up 
to 1990 have been summarized, for example, by van Kam- 
pen Q . In general such efforts have been regarded as not 
very successful 10]. 

In spite of the difficulty, there have been continuous ef- 
forts on the construction of cost function and related top- 
ics. Elegant results have been obtained in several direc- 
tions. Tanase-Nicola and Kurchan [ll| have considered 
explicitly the saddle points of gradient systems. They 
started from the existence of potential or cost function 
to avoid the most difficult problem of the irreversibil- 
ity. The gain is that they can now obtain a powerful 
computational method to count the saddle points and to 
compute the escape rate. They also provide an extensive 
list of related literature. 

The study on the mismatch of the fixed points of the 
drift force and the extremals of the steady state distribu- 
tion has been reviewed by Lindner et al. |l2j . Rich phe- 
nomena have been observed, but the mismatch has been 
treated as "experimental" result. There is no mathemat- 
ical/theoretical explanation on why it should happen. 

In another survey the useful and constructive role 
played by the noise has been demonstrated by exam- 
ples ^3|- ^ ls argued that the noise is essential to es- 
tablish the functions of dynamical systems. Again, the 
mismatch problem is encountered and the constructed 
potential function is often regarded as approximation. 

From a different perspective, there has been an effort to 
provide a solid foundation for non-equilibrium processes 
based on the chaotic hypothesis |l4|. The chaotic hy- 
pothesis presumes that the system is sufficiently chaotic 
that variation of parameters of the system leads to a 
unique parameter-dependent steady state, even though 
the Gibbs entropy change is not path independent |l5j|. 
Under this hypothesis an interesting and important fluc- 
tuation theorem has been obtained, which further sug- 
gests the existence of the Boltzmann-like steady state 
distribution function. Hence, a cost function very likely 
exists under this situation. A difficulty with this ap- 
proach is that extremely few practical physical systems 
have been shown to satisfy the chaotic hypothesis. 

Because the metastability is such an important phe- 
nomenon and because of the difficulty encountered in the 
construction of cost function, efforts have been made to 
go around the cost function problem when computing the 
life time of a metastable state. The effort results in the 
so-called Machlup-Onsager functional method, summa- 
rized by Freidlin and Wentzell llr|. T his approach has 
been actively pursued recently |l7l.ll8| . 

In this paper we give a careful discussion of the basis of 
our method of constructing a cost function in the linear 
region close to any fixed point, whether stable or unsta- 
ble. We show that this approach gives an unambiguous 
prescription under wide conditions. The only case we 
have found which does not give an unambiguous expres- 
sion for the cost function has a subspace within which 
the noise does not act, and out of which the force does 
not carry the state. 



In section [H] we give a general discussion of a linear 
system with noise. In section lllll we show how a cost 
function can be constructed by a decomposition of the 
force matrix into two factors, one of which is a symmet- 
ric cost function matrix, and that this gives a probability 
density of exponential form, with the exponent propor- 
tional to the cost function. The general proof of this 
result is obtained in the Appendix, where we exploit the 
Jordan transformation for matrices with an incomplete 
set of eigenvectors. In section we discuss the singu- 
larities of the transformation, and identify two types of 
singularities, one of which corresponds to a flat subspace 
of the cost function, while the other corresponds to the 
possibility that the dynamics separates the system into 
two or more disjoint subspaces. In section we discuss 
solutions other than the Boltzmann-like solution of the 
equation for a stationary distribution, and argue that 
such solutions may be significant in any attempt to ex- 
tend this solution to the nonlinear regime. In section IVT1 
we discuss the significance of this decomposition of the 
force matrix, and its relation to the principle of detailed 
balance. 



II. 



STATING THE PROBLEM 



Many processes in natural sciences can be modeled 
quantitatively. One particularly important class of such 
modeling is that described by first order differential equa- 
tions pj| , supplemented by stochastic terms Q ■ We start 
with the nonlinear dynamic equation 



x = f{x) + C(x,t) 



(1) 



which gives the stochastic evolution of the state repre- 
sented by the real d-dimensional vector, x T = (x\, Xd)- 
Here the superscript r denotes the transposed vector. 
The force vector is the d-component f(x), and this gives 
the deterministic time evolution of the system. For sim- 
plicity we take the noise term t) to be Gaussian 
white noise, with zero mean, (C(t)) = 0, and variance 
{C(x,t)C{x,t')) = 2D(x)S(t - t'). The angular brackets 
denote the average over noise distribution, and 5(t) is the 
Dirac delta function. In this work we assume D(x) to be 
independent of x. The probability distribution function 
p(x, t) then satisfies the Fokker-Planck equation 



dp(x,t) 
dt 



d 

dxi 



~Mx)+D^ 



p(x,t) . (2) 



In the neighborhood of a fixed point, which we take to 
be at the origin, the force can be replaced by its linear 
approximation 



fi{x) -^ijXj 



(3) 



It was noticed by Ao [|j that, in the linear region near a 
stable fixed point, Eqs. Q and can be decomposed 
in the form 



(S + A)x = -Ux + £(t) 



(4) 
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where the symmetric matrix S is semi-positive definite 
and the matrix A is antisymmetric. The noise function £ 
has variance given by t)^ T (x, t')) = 2SS(t — t'). In 
this paper we adopt an equivalent, but simpler, approach 
of factorizing the force matrix as 

F= -{D + Q)U= -{S + A)~ l U , (5) 

where D is the symmetric diffusion matrix, Q is an an- 
tisymmetric matrix which can be determined, and U is 
the symmetric cost function matrix, which was called a 
potential matrix in ref. This breaks the force Fx 
into two components, F^x — —DUx, which generates 
a motion towards the origin if U is positive definite, and 
F ( c ) x = -QUx, which gives a motion on the manifold of 
constant U(x) = (l/2)x T Ux. The quadratic form li(x) is 
the cost function. 

It is immediately obvious that if the vector / is re- 
placed by F^ d 'x = —DUx in Eq. J2J), a Boltzmann-like 
stationary distribution of the form: 

p{x) oc exp{— U{x)} (6) 

satisfies the equation, since the current density 

ji(x) = [fi(x) - Dijd/dXj] p (7) 

vanishes. Since 

F {c) 

x generates a current density 

j^(x) = -QUxp(x) , (8) 

which is divergence free, and conserves IA and p, the com- 
bination fix) = [F^+F^a; also conserves the distribu- 
tion given by Eq. JSJ, so that this is a stationary solution 
ofEq. ©. 

This decomposition of the force matrix allows an ex- 
plicit time-independent solution of the Fokker-Planck 
equation to be written down. In ref. [S] the solution of 
the equation for Q was obtained by a power series ex- 
pansion, without any discussion of the convergence of 
this series. In the next section we show that there is a 
unique solution for the equation for Q under rather wide 
conditions. It is not even required that the fixed point 
of F be stable, although if it is not the stationary solu- 
tion given by Eq. (jSJ) is unbounded, and could only, at 
best, give a useful solution in a neighborhood of the fixed 
points with boundary conditions that do not perturb this 
solution too strongly. 

It is worth noting that, under the coordinate transfor- 
mation x — > y = M x, the mapping F that relates x to 
x transforms as 

F -> M^FM , (9) 

while the symmetric matrices U and D, which represent 
quadratic forms, transform as 

U^M T UM, D —> M^D (M T )~ 1 , (10) 

and Q transforms in the same way as D. These transfor- 
mations preserve the separation of D + Q into symmetric 
and antisymmetric parts. 



Starting from the work of Onsager Q , there has been 
an extensive literature on dynamical behavior near a sta- 
ble fixed point 20]. The new construction clearly works 
for this important situation, and has indeed offered a 
new angle. However, it is not sufficient to generalize the 
decomposition of the force in Eq. (J5J to the nonlinear 
regime by an equation of the form / = — [D + Q)gra,dU, 
since the Boltzmann form, Eq. © may not satisfy the 
Fokker-Planck equation (j2J if the antisymmetric matrix 
Q is space dependent. The generalization to nonlinear 
systems needs further study. In the next three sections we 
establish the decomposition firmly in the linear regime, 
and investigate its limitations and implications. 

III. DECOMPOSITION OF THE FORCE 
MATRIX 

In this section we develop a general method for making 
the decomposition of the force matrix given in Eq. JSJ. 
Since this equation, together with the symmetry of U, D, 
and antisymmetry of Q, leads to 

U = -(D + Q)- 1 F = -F T iD-Q)- 1 , (11) 

the equation to determine Q is 

FQ + QF T = FD — DF T . (12) 

This is a system of did — l)/2 linear equations to deter- 
mine the same number of independent components of Q, 
so it has a unique solution unless the set of equations is 
singular. Inversion of the matrix D + Q then gives the 
matrix S + A of Eq. |QJ. 

Our method of solution is best illustrated by consid- 
ering the case that F is real symmetric or has distinct 
eigenvalues, so that it can be diagonalized in terms of its 
left and right eigenvectors. Equation l|12|l then takes the 
form 

(A Q + \p)Q a (3 = (A a — ^p)F) a p , (13) 

where the A Q are the eigenvalues of F, and the tilde de- 
notes this representation in terms of eigenvectors. This 
gives an immediate solution for Q provided that no pair 
of the eigenvalues of F adds up to zero. The eigenvalues 
can only add to zero when the fixed point is unstable, 
and this is discussed in section llVl 

For completeness we must consider the general case 
with degenerate eigenvalues for asymmetric F, in which 
case there may not be a complete set of eigenvectors. 
This case is dealt with in the Appendix, using the Jordan 
representation of a nonsymmetric matrix. 

IV. SINGULARITIES OF THE 
TRANSFORMATION 

There are two places in our argument where the trans- 
formation from the force matrix F to the symmetric cost 
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function matrix U might be singular. Equation l|12[l for 
the antisymmetric matrix Q can be solved, and we have 
an explicit solution in section ITTll unless the determinant 
of the coefficients in d(d — l)/2 inhomogeneous equations 
is zero. The second possibility is that the matrix D + Q 
whose inverse appears in Eq. might have zero deter- 
minant. 

In section IIIII we showed that the conditions for the 
equation for Q to be singular are that two of the eigen- 
values of F sum to zero, or, as can be seen from Eq. (IA8() 
in the Appendix, when the null space of F 2 has two or 
more dimensions. Neither of these cases arise for a stable 
fixed point. There are two distinct cases of X a + Xp = 0, 
according to whether the two eigenvalues are real eigen- 
values of opposite sign, or whether they form a complex 
conjugate pair. We study the behavior of the eigenvalues 
and eigenvectors of U in these two cases, assuming that 
the two eigenvalues of F are nondegenerate, and that F 
has no zero eigenvalue. 

Instead of studying the eigenvectors of U directly, we 
study the eigenvalues and eigenvectors of 

U- 1 = -F-\D + Q) = -RA-\D + Q)R T , (14) 

which has the same eigenvectors but reciprocal eigenval- 
ues. Here R is the matrix whose columns are the right 
eigenvectors of F, and A is the diagonal matrix with the 
eigenvalues of F as its diagonal elements. The general- 
ization of these definitions of R and A to the case where 
the eigenvalues of F are not complete is given in the Ap- 
pendix in Eqs. 1|A3(1 . (|A5|) and ijAlOfl . For a pair of eigen- 
values with X a + Xp 0, with no other sums of two eigen- 
values small and no other small individual eigenvalues, 
the only large terms in the matrix U^ 1 = —A~ 1 (D + Q) 
are, according to Eq. (JT^J, 



This gives the two small real eigenvalues of U as 



2D 



a/3 



X a + A/} 



(15) 



When all other matrix elements are neglected, this ap- 
proximation, combined with Eq. I|14(l . gives 



X a + A/3 



2D, 



a/3 



R ja R ]fi ± /X! R i<* X! R . 
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and the corresponding eigenvectors as 



22 « ±. 
ap l 



2j Rja 



(19) 



(20) 



For the case of a pair of real eigenvalues of opposite 
signs we can see that, as the sign of X a + Xp is changed 
by a change in the parameters of F, D, the stable and 
unstable manifolds of U change places with one another. 
The stable and unstable manifolds of U bisect the stable 
and unstable manifolds of F in the original representa- 
tion, as is obvious if one normalizes the real eigenvectors 



of U by R la = Ria/ JJ2j R ja and Riff similarly, giving 

eigenvectors, Ri a ± Rip. In the limit X a = Xp, U is flat 
in this two-dimensional subspace. 

For a complex conjugate pair of eigenvalues with 
A Q + A* w the behavior is a little different. Using 



Rjp, we find the term inside the 



the property Rj a 
bracket in Eq. I|19f) is always positive. The two eigenval- 
ues of U then have the same sign, so U is either stable or 
unstable, depending on the sign of X a + Xp, in this two di- 
mensional subspace. As the real part of X a changes sign, 
a two-dimensional stable manifold becomes unstable, or 
vice versa. 

It can be seen from Eq. (|A7() in the Appendix that 
where one of the eigenvalues satisfying A Q + A^ = cor- 
responds to a higher dimensional subspace there may be 
higher order zeros of the eigenvalues of U. 

If the matrix D is positive definite there is no possibil- 
ity that D + Q could be singular. If u is a vector in the 
null space of D + Q, we have 



= u T (D + Q)u = u T Du, 



(21) 



2D a p 
X a + Xp 



(S la Rip + 8-ipRia) ■ (16) 



This equation, combined with the relation RL = /, 
shows that the approximate eigenvector corresponding 
to a large eigenvalue w _1 can be written as 

a a R ia + o.pR l p = y^^agRja + apRjp) R^L^j , (17) 



provided the amplitudes and eigenvalues satisfy the equa- 
tion 



since the antisymmetry of Q makes its expectation value 
vanish. Therefore D ± Q cannot be singular when D is 
positive definite. 

However, we do not usually want to specify that the 
noise acts on all coordinates. Typically, when two of the 
coordinates are the position and momentum of a particle, 
people will take the noise to change the momentum but 
not the position of the particle. However, Eq. (|21|l shows 
that for non- negative definite D, vectors in the null space 
of D + Q are in the intersection of the null spaces of D 
and Q. Equation (|12fl then shows that, for such a vector 
u in the null space of D and Q, 



'30, 



W 1 ap - - x^+Xfi ( 53j R 2 ja + a (i Hj R 3aRj0 



= F(D - Q)u =(D + Q)F T u , 



(22) 



and so u is only in this null space if F T u, or any power 
(18) of F T acting on u, is still in the null space. 
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This condition is in agreement with what one should 
expect. The noise does not have to act directly on all 
coordinates, but, if there is a subspace in which there 
is no noise, and which is left invariant by the motion, 
there can be no equilibration within that subspace except 
collapse towards a stable fixed point. 

V. OTHER STATIONARY SOLUTIONS 

Although the Boltzmann-like form given in Eq. © 
gives a stationary solution of the Fokker-Planck equa- 
tion, it is only the unique solution under certain rather 
restrictive boundary conditions. One can see clearly why 
this might be an issue by considering the one-dimensional 
form of the equation near a stable fixed point, which can 
be written as 

S + i^» = °- < 23 » 

In addition to the Boltzmann-like solution p(°~) — 
(27r) -1 / 2 cxp(— x 2 /2), this has a current-carrying solution 
of the form 

pW(x) oce-* 2 / 2 [ X e x ' 2 / 2 dx' . (24) 
Jo 

This expression is proportional to 1/x for large values of 
x, so, if the linear approximation to the equation is valid 
up to fairly large values of x, the coefficient of such a term 
must be exponentially small to prevent the probability 
density given by p^ + p^ from being negative. 

In d-dimensional systems there are similar solutions 
falling off like l/|x| d for large |x|. For such solutions of 
Eq. J21 the current at large distances from the origin is 
primarily driven by the linear force Fx, and the diffusive 
motion is a small correction, so, while the density falls 
off like |x|~ d , the conserved current falls off like |x|~ d+1 . 
Again, current conservation shows that this contribution 
to the density must be positive and negative in different 
parts of space, so that its coefficient at the origin must be 
exponentially small, with an exponent that depends on 
the size of the region in which the linear approximation 
is valid. 

Near a maximum of the cost function, where p^ oc 
exp(x 2 /2), the one-dimensional current-carrying solution 
has the form 

p«(x) oce* 2 / 2 I* e- x ' 2 ' 2 dx' . (25) 
Jo 

This grows at large distances in the same way as p^ ', 
with a change of sign at the origin. In d dimensions a 
saddle point can sustain a relatively large current across 
it, because there are current-carrying states for which p 
is of the same order of magnitude as p^ . 

To get such a current across a stationary point in a lin- 
ear system, it is necessary to impose some external cur- 
rent sources and sinks. However, if we want to describe a 



nonlinear force field in terms of its approximately linear 
behavior in the neighborhood of its fixed points, adjacent 
neighborhoods can generate external current sources and 
sinks for one another, so we should not be surprised to 
find such currents if we linearize in a local region. These 
"external" currents will not only produce flows at the 
boundaries, but will shift the flow lines away from the 
surfaces of constant cost function U shown in Eq. JSJ. 

In the neighborhood of the minimum of the cost func- 
tion, such current-carrying solutions resembling the one- 
dimensional Eq. I|24l) shift the maximum of the den- 
sity away from fixed point, since the gradient of p^ is 
nonzero. Since, as we remarked in connection with this 
equation, the amplitude of such a term must fall off ex- 
ponentially with the size of the region of linearization, 
in order to prevent negative densities, it should not be 
possible to obtain such a term by a conventional per- 
turbation theory in the neighborhood of the stable fixed 
point. Our numerical exploration of nonlinear systems 
of this sort suggests that these current-carrying solutions 
are significant, since the maximum of the density is dis- 
placed from the zero of the force. One possibility is to 
introduce such current-carrying states, in addition to a 
p^ determined by the cost function, in order to make 
this method applicable to nonlinear systems. 



VI. DISCUSSION 

The main result of this work is to show that, for a sys- 
tem with a deterministic motion controlled by a linear 
force, and a diffusive motion driven by constant white 
noise, the force matrix F can be decomposed into two 
parts, —DU and —QU. Provided the fixed point of F 
is stable, the first of these components leads to a steady 
state distribution of the Boltzmann form, exp(— ti(x)), 
with no probability current, the usual form of an equi- 
librium distribution when detailed balance holds. The 
second part gives a flow on the surface of constant U. 

The cost function matrix U can be diagonalized by an 
orthogonal transformation, and if it is positive definite, 
it can be transformed to the identity matrix by choosing 
a new scale for the variables. In this representation the 
dissipative part of the force matrix is F^ d > = —D, which 
is suggestive of Einstein's relation between diffusion and 
dissipation , or of the fluctuation-dissipation theorem 

m 

When the cyclic motion induced by F^ is included the 
relation between the eigenvalues of F and D becomes 
more complicated. If the motion in a two-dimensional 
subspace is dominated by a fast cyclic motion within the 
subspace, there will be a complex conjugate pair of eigen- 
values of F, so that they have a common relaxation rate 
given by the real part of the eigenvalues. These possibil- 
ities require much more detailed work than we have yet 
given them. 

This is actually not so different from the situation usu- 
ally encountered in statistical mechanics, at least if a 
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classical rather than a quantum description is used. For 
a damped harmonic oscillator there is a cyclic motion 
in phase space, as well as the thermal noise and viscous 
damping acting on the momentum coordinate, while for 
black-body radiation in a cavity there is cyclic motion 
between electric and magnetic fields, in addition to the 
resistive damping and noise from the walls of the cavity. 

What is remarkable is not that the steady state den- 
sity can be written as the exponential of a cost function, 
since if there is a steady state we could always define the 
cost function as minus the logarithm of the steady state 
density. We find it remarkable that for a linear stochas- 
tic system of this sort it is generally true that the force 
can be decomposed into two parts, one of which gives de- 
tailed balance in its strictest sense, while the other gives 
a cyclic motion on the surfaces of constant cost function. 

The question of whether this technique can be ex- 
tended to nonlinear systems is an important one, but 
requires careful investigation. We have done preliminary 
work based on perturbative inclusion of nonlinear terms, 
and made comparison with numerical calculations. It is 
clear that Eqs. (0J and (J5J need some modification for 
nonlinear effects, and also that if we try to solve the 
problem by matching limited regions in which linearity 
holds approximately, the matching may introduce solu- 
tions of the linearized Fokkcr-Planck equation other than 
the Boltzmann-like p(°K Different regions will have to 
serve as sources and sinks for their adjacent regions. 
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APPENDIX A: APPENDIXES 

For cases in which the force matrix F has degener- 
ate eigenvalues and is not symmetric, so that it may not 
have a complete set of eigenvectors, we use the Jordan 
transformation psj of a general real matrix. The Jor- 
dan transformation uses a complete set of independent 
column vectors v a with the property 



Fl! a = X a V a + [i a ~!V C 



(Al) 



where /i Q is zero if X a =/= A Q+ i and is either unity or 
zero for A Q = A a +i. The set of row vectors u a with the 
orthonormality property u a v@ — S a p then satisfies the 
equation 



U a F = X a u a + p a u aA 



(A2) 



Let us define matrices R and L as Ri a = vf and L a i = 
uf . Since the vectors u a are orthonormal to the set v a , 
we have LR = I, so that these matrices are inverse to 
one another. Ri a , L a i are real for real X a . For complex 
X a , its complex conjugate is also an eigenvalue, Xp = A* , 
since F is a real matrix. In this case R a i = R*p^ Li a = 
L*o. The Jordan transformation is then given by 



LFR = A 



(A3) 



The matrix A is now block diagonal, where each nonzero 
block has identical diagonal elements, which are degen- 
erate eigenvalues, and unity in each place immediately 
above the diagonal. 

With this result, Eq. I|12l) can be rewritten as 



where 



AQ + QA T = AD - DA T , 



= LQL T , D = LDV 



(A4) 



(A5) 



The matrix Q remains antisymmetric and D remains 
symmetric. In this representation Eq. (| A4K takes on the 
form 

a, 0+1 

= (A Q - Xfj)D a [} + ii a D a+ i.j3 — ppD a ^ + i . (A6) 
It has the solution, for a < f3, 

~ _ Aa-A/3 ~ 
Wa/3 — ~ —^—J-^af3 



Aq + Xf 



-2A, 



^ (x a + x^ Da+ ^ 



fj,>l 



-2X a Y, 



(-iy 



U>1 



(Aq + X/3) 



V+ 



-2E 



^! (Aq + XpY + U+l 



(A7) 



where the sums go over all values of /i, v for which the 
indices lie within the same block of the block-diagonal 
matrix A. For the case A Q = Xp this reduces to 



sr n i \- (-1)" n 



(n + v-i)i {-iy +v 



^ [i\v\ (2Aq)^+ 



For the case in which all the fi a are zero, Eq. (| AT|) is 
equivalent to Eq. (|13fl . 

In this Jordan representation of force matrix F, the 
state variable is transformed by 



y = Lx , 



(A9) 
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and the corresponding transformation of U is given by which is form-invariant with Eq. under the Jordan 

transformation. 

U = R T UR . (A10) 
Equation then becomes 

(D + Q)U = -A, (All) 
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